m 
o 
o 

(N 



DUKE-TH-03-237, INLO-PUB-02/03 

Chiral Limit of Strongly Coupled Lattice Gauge Theories 

David H. Adams"'^ and Shailesh Chandrasekharan" 
" Department of Physics, Box 90305, 
Duke University, Durham, 
North Carolina 21108, USA 
and 

^ Instituut-Lorentz for Theoretical Physics, 
Universiteit Leiden, Nield Bohrweg 2, 



; 2300 RA Leiden, The Netherlands. 



Abstract 

We construct a new and efficient cluster algoritlim for updating strongly coupled U{N) lattice 
gauge theories with staggered fermions in the chiral limit. The algorithm uses the constrained 
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. monomer-dimer representation of the theory and should also be of interest to researchers working 

, on other models with similar constraints. Using the new algorithm we address questions related 

m ' 

' to the chiral limit of strongly coupled U (N) gauge theories beyond the mean field approximation. 

. We show that the infinite volume chiral condensate is non-zero in three and four dimensions. 

Qh! However, on a square lattice of size L we find '0V'(O)) ~ L'^"'^ for large L where 

' r/ = 0.420(3)/A^ + 0.078(4)/A^^. These results differ from an earlier conclusion obtained using a 

, i ■ different algorithm. Here we argue that the earlier calculations were misleading due to uncontrolled 

^ , autocorrelation times encountered by the previous algorithm. 
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I. INTRODUCTION 



One of the challenges in the study of strong interactions is to compute physical quantities 
with in the framework of QCD with controlled errors. Although lattice QCD can in principle 
accomplish this task, most known algorithms encounter problems related to critical slowing 
down as the quark masses become small It is impossible to compute quantities reliably 
for realistic "up" and "down" quark masses with current algorithms. At present, one usually 
computes quantities at an unphysically large quark mass and uses chiral extrapolations to 
obtain the final answer P| . For sufficiently light quarks this is a reliable technique. However, 
for the quark masses that are currently used, such extrapolations are questionable. Further- 
more, most calculations with light quarks are obtained in the quenched approximation where 
the internal dynamics of fermions are ignored. There also exist some interesting quantities 
which cannot be computed when the quarks are heavy. For example, in the real world the 
rho meson can decay into two pions. Unfortunately, for most values of the quark masses 
that are currently used, this decay is forbidden. If lattice simulations were possible with 
sufficiently light quarks, the resonant nature of the rho meson can be studied with lattice 
techniques j3|. 

There are two main problems in constructing efficient algorithms for lattice QCD. The 
first is that the fundamental variables, namely the quark and gluon fields, are unphysical 
since they are gauge dependent. Secondly, the quarks are fermions and due to the Pauli 
principle introduce sign problems which are typically very difficult to solve. Fortunately, 
with a good lattice fermion formulation and at zero baryon chemical potential, the second 
problem can be avoided. The quarks can be integrated out and their effects can be encoded 
in terms of a determinant which is a positive non-local function of the gauge fields. However, 
one is still left with the problem of updating the unphysical gauge fields. Previous experience 
suggests that in order to find an efficient algorithm it is often useful to understand and 
isolate the physical modes of the theory. Unfortunately, for gauge theories, especially with 
the non-local fermion determinant this problem looks formidable. 

Interestingly, there is a limit of lattice QCD with staggered fermions where most of the 
above mentioned complications can be eliminated. This is the so called strong coupling limit 
^. Although this limit has the worst lattice artifacts and could describe the wrong phase, 
the resulting theory is still an interesting toy model for QCD. It contains the physics of 
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confinement and chiral symmetry breaking. In the chiral limit one finds massless Goldstone 
bosons in addition to other massive mesonic and baryonic excitations. The strong coupling 
theory can be solved analytically only in the large N and large d hmits. One needs a 
numerical approach to study it without these approximations. 

The luxury of the strong coupling limit is that it is possible to integrate out the gauge 
fields analytically and write the problem entirely in terms of gauge invariant objects. It was 
pointed out in 7| that with U{N) gauge fields the strong coupling model with staggered 
fermions is equivalent to a system of monomers and dimers with positive Boltzmann weights. 
Later a monomer-dimer-polymer representation was discovered for SU{3) gauge fields Is 
it possible to use these simplified representations that arise at strong couplings to construct 
an efficient algorithm in the chiral limit? In simple local algorithms were proposed. 

However, as far as we know, all these algorithms break down in the chiral limit. In fact there 
is evidence that the proposed algorithms may also have other problems away from the chiral 
limit y. Recently, a cluster based algorithm was proposed to update the monomer-dimer 
system which works very well in the chiral limit on small lattices Q, However, even 
this algorithm becomes inefficient for large lattice sizes due to long auto correlation times. 

In this paper we propose a new and efficient algorithm to study the chiral limit of strongly 
coupled U (N) gauge theories with staggered fermions in any dimension. It should be possible 
to extend it to SU{N) gauge theories with minor modifications. We then use this algorithm 
to study questions related to chiral symmetry in two, three and four dimensions. Our study 
shows that this algorithm has the potential to address many unresolved questions about the 
chiral limit of gauge theories at least in the strong coupling limit. Our paper is organized as 
follows. In section 2 we review the monomer-dimer representation of strong coupling U{N) 
gauge theories with staggered fermions and discuss the consequences of chiral symmetry in 
this language. We then show how a finite size scaling formula for the chiral susceptibility 
can be used to compute the chiral condensate from finite volume lattice calculations in the 
chiral limit. In section 3 we construct a new type of cluster algorithm, which we call the 
"directed path" cluster algorithm, to update the monomer-dimer model. We show that 
it obeys detailed balance and is ergodic. Section 4 contains explicit expressions that can 
be used to compute observables like the condensate and the susceptibility with the new 
algorithm. In section 5 we test the algorithm by comparing exact results on small lattices 
with the results obtained using the algorithm. We also present an exact large L result on 
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the 2 X L lattice and compare with the result from the algorithm. In section 6 we discuss 
;he performance of the new algorithm and compare it with an earlier algorithm proposed in 
10 . h|. We argue why the results obtained earlier were incorrect at large volumes. In this 
context we also briefly study the autocorrelations of the new algorithm. In section 7, using 
our new algorithm we study the issue of chiral symmetry breaking in two, three and four 
dimensions for different values of A^. Section 8 contains our conclusions. 



II. MONOMER-DIMER REPRESENTATION 

Let us review the monomer-dimer representation of the strongly coupled U{N) lattice 
gauge theory with staggered fermions. The Euclidean space action of the model we consider 
is given by 

S = —"^i^fj^^x) %jj{x)U^{x)^p{x + jl) — 'ijj{x + p.)Uj^{x)^{x) — m'^'ijj{x)ip{x) (1) 

X,fJ, X 

where x = {xi,X2, Xd) labels the sites on a d-dimensional periodic, hyper-cubic lattice and 
H = 1,2, ..,d labels the various directions. For concreteness we assume x^ G 0, 1, 2, — 1 
such that is the length of the hyper-cubical box in the fi direction. We will use V = 
L1L2 ■ -Ld to denote the volume of the lattice. The site next to x in the positive /i direction 
is labeled x + /t. The link variables connecting the sites x and x + fi, represented by U^{x), 
are N x N unitary matrices. ip{x) is an component column vector and is an 

component row vector. Both these vectors are made with Grassmann variables and represent 
the staggered fermion fields at the site x. We will assume that the gauge links satisfy periodic 
boundary conditions while the fermion fields satisfy either periodic or anti-periodic boundary 
conditions. The factors rjnix) are the well known staggered fermion phase factors. However, 
we will choose them to be rii{x) = t and rifj_{x) = exp[z7r(a;i + .. + a;^_i)],yU = 2,3, ...d, 
where t is a real-valued parameter incorporating the effects of temperature. By working on 
„t„c lattices with L, « = 2,3 . and allowing < to va.y continuously, we can 

study finite temperature phase transitions [1^. When t = 1, the ri^{t) turn into the usual 
phase factors. The parameter m controls the fermion mass. Our definition of the action, 
(eq.(^), differs from the conventional definition since we have dropped a factor of half in 
front of the kinetic (hopping) term. This only changes the normalization of the fermion 
fields and the definition of the fermion mass m up to a factor of 2, but helps in avoiding 
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extra powers of two in many expressions we write below. 
The partition function of the model is given by 



n [dU^{x)] n dij{x)] exp(-5), 



(2) 



where [dU] is the Haar measure on the U{N) group, [dilj{x)d'ilj{x)] represent Grassmann 
integration. The model with t = 1 was considered earlier in IJ], where it was shown that 
the integrals in eq.(jD) can be performed analytically and the partition function can be 
rewritten as a sum over positive definite Boltzmann weights associated to monomer-dimer 
configurations. Let us see this explicitly by setting t = 1. First note that the integral over 
the gauge fields can be done one link at a time in the background of the Grassmann fields. 
In LZll it was shown that 



[dU^{x)] exp ^{x)U^{x)ip{x + fi) — ip{x + jl){U^{x))'^ip{x) 



^ {N-b)\ 

b=Q 



N\b\ 



■ip{x)4'{x)ip{x + jl)ip{x + fl) 



and 



[dil){x) dil){x)\ exp{mip{x)ilj{x)) i/j{x)i/j{x 

Using these relations we can rewrite the partition function as 

{N-b,ix)) 



(N-n) 



-m 



nl 



E n 

InM 



b^ix)\N\ 



TT ^' r 

\ nJ 



(3) 



(4) 



(5) 



0, 1,...,A^ is the number of monomers located at the site x and bfj,{x) 



where Ux 

0,1,2, N is the number of dimers located on the bond connecting the sites x and x + fi. 
The configurations [n, b] denote the sets of monomer values n = {Ux} on all sites and bond 
values b = {bf^{x)} on all links which satisfy the constraint 



nx + E 



b^^{x) 



N 



(6) 



on each site x. We use the definition b-^{x) = b^{x — fi). In other words the total number of 
monomers and dimers associated to a site must always be A^. The only effect of a general t 
is that every dimer along the /i = 1 link is weighted by an extra factor of t^. For illustration, 
an = 3 monomer-dimer configuration is shown in figure 
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FIG. 1: An example of N = 3 monomer- dimer configuration on a A x 4 lattice. 



When m = the action given in eq. is invariant under U{1) chiral transformations, 



'(/'(x) e^^il){x), ijj{x) —> ip{x)e^^, for x even 



■i/j{x) 



^{x), ijj{x) — * ip{x)e , for x odd. 



(7) 



A site X is defined to be even (odd) if xi +X2 + ... + is even (odd). To study the properties 
of the vacuum under chiral transformations one usually defines the chiral condensate, 



(8) 



In a finite volume the vacuum is chirally symmetric in the chiral limit. A consequence of 
this symmetry is that 

lim {ijjilj) = 0. (9) 



m— >0 



In the monomer-dimer model 



1 

mV 



(10) 
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where {rix) is the average number of monomers on the site x. In the chiral hmit at finite 
volumes this average goes to zero as m^, which again means the chiral condensate vanishes. 
Mean field calculations on the other hand suggest that the strong coupling vacuum breaks 
the chiral symmetry spontaneously [s,!^. This means 

lim hm (^^) ^ (11) 

We will call this non-vanishing quantity the infinite volume chiral condensate. It would be 
interesting to calculate this condensate using the new algorithm. However, using eq. ()ll|) to 
extract it is difficult since one needs to add a small fermion mass term to the action and 
then perform a careful scaling analysis of the chiral condensate with respect to both the 
volume and the mass. Since our algorithm can work efficiently at m = in a finite volume 
we will instead calculate the chiral susceptibility, defined by 

It is easy to show that 

^-"'(<«)^^ (13) 
Usually, the chiral susceptibility is defined as the first derivative of the chiral condensate with 
respect to the mass. However, our definition includes the disconnected term proportional to 
the square of the infinite volume chiral condensate. With this definition, in a cubical box 
\^ = L'^, we can argue that 

const. L'^ in a phase with broken chiral symmetry 
X = ■( const. L^"'' at a chirally symmetric critical point (14) 
const. in the chirally symmetric massive phase 

in the large volume limit. In the phase with broken chiral symmetry the coefficient of L*^ 
is the square of the infinite volume chiral condensate which can be extracted by studying 
the large volume behavior of x and fitting the data to the above form. With conventional 
algorithms x is very difficult to compute since it is a noisy observable. Our new approach on 
the other hand allows us to measure it very accurately even in large volumes. Although we 
have written eq.(|14p for a general c?, the Mermin-Wagner-Coleman theorem 3 forbids 
a phase with broken chiral symmetry in two dimensions. In that case only the critical or 
the massive phases are possible. 



III. DIRECTED PATH ALGORITHM 



The constraints imposed by eq. (jHl) make it difficult to design algorithms for the monomer- 
dimer systems near the chiral limit. When m = 0, configurations with monomers have 
vanishing weight in the partition function and local algorithms find it difficult to update 
the remaining constrained dimer configurations efficiently. Cluster algorithms on the other 
hand can deal with constraints very efficiently. An analogue of this problem in a well-known 
setting is the following: Consider for example a quantum spin-half system where the number 
of "up" spins and "down" spins are individually conserved. A typical configuration of such 
a system is represented by a world-line of spins. While local algorithms find it extremely 

n 

difficult to update such configurations, the loop cluster algorithm is very efficient J^]. In 
lol we discovered that the monomer-dimer model can also be rewritten in terms of 
loop variables and found them to be convenient tools to update the system when m = 0. 
Unfortunately, the Metropolis algorithm that was designed to update the loops turns out 
to be inefficient. Certain large clusters can only be updated with small probabilities and 
the algorithm develops long autocorrelation times as the lattice size grows. This problem 
is similar to the problem encountered in an anti-ferromagnetic quantum spin-half model in 
the presence of a uniform magnetic field. The loop cluster algorithm which works efficiently 
in updating the model in the absence of magnetic fields, becomes exponentially inefficient 
in their presence. Again, some large clusters get frozen and cannot be updated. 

Recently, a new algorithm called the "directed loop" algorithm was proposed for the 
antiferromagnetic model in a magnetic field This algorithm is extremely efficient even 
for large magnetic fields. Here we extend this algorithm to study strong coupling gauge 
theories. In this section we construct a "directed path cluster algorithm" to update the 
monomer-dimer configurations. Our construction is such that the algorithm does not change 
the number of monomers in a given configuration. Hence, it is only ergodic in the chiral 
limit where the number of monomers is strictly zero and does not change. When m 7^ it 
is necessary to supplement it with an update that changes the number of monomers. For 
this we can use any local algorithm. We will briefiy review one such algorithm at the end 
of this section for completeness. 

The basic idea behind our algorithm is to create a monomer at a site. However, in order 
to do this while satisfying the constraints it is necessary to remove a dimer from one of 
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the links connected to the site and create another monomer on the other end of the bond. 
This neighboring monomer is then moved along a directed path, while satisfying detailed 
balance at each step until it encounters a partner such that both monomers can be removed 
by creating a dimer. Thus, at the end of a directed path update, the number of monomers 
remains fixed. It is interesting to note that at every stage, between the start and the end 
of the path, we sample configurations that are relevant in the computation of two point 
correlation function of monomers. Hence, these intermediate configurations can be used 
to measure certain observables. We will discuss this feature of the algorithm in the next 
section. 

In order to explain the algorithm better and show that it obeys detailed balance it is 
useful to develop some notation. We begin by dividing the sites of the lattice into active 
and passive sites. If the first site we pick during the update is even then all even sites are 
defined to be active and all odd sites become passive. On the other hand if the first site 
happens to be odd, then all odd sites become active and even sites become passive. Each 
active (passive) site is associated with an active (passive) block which includes the site with 
the information n about number of monomers on it and the 2d bonds connected to the site 
with the information 61, ...6;^, 6_d of their dimer content. Due to the constraint eq.ijHl) 
we must have 

n + bi + b^i + +... + bd + b_d = N. (15) 

We will represent the block symbolically by (n, ...,bd,b_d)- Its Boltzmann weight is 

defined to be 

~ n\ l[ b,ib.,\imy 

if it is an active block and 

W^p™ = ^ (t2)(^^+^-0 (17) 
nl 

if it is a passive block. We have dropped factors of the mass because our objective is not 
to change the total number of monomers in the configuration. It is easy to verify that the 
product of the weights of all the active and passive blocks in a configuration is equal to the 
Boltzmann weight of the configuration up to some power of the mass. 

A directed path update begins by entering a site at random. By definition that site 
belongs to an active block and the path enters the block through the site. Given an in- 
coming path the algorithm works on finding an outgoing path with a probability that satisfies 
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detailed balance. The out-going path can either be one of the 2d bonds or the starting site 
itself. As we will see shortly, the correct probability to choose the outgoing path to be the 
starting site is 

Pss = ^ (18) 
and the probability to choose the bond in the /U = 1, —1, 2, —2, .., d, —d direction is 

Ps, = |. (19) 

If the out-going path is chosen to be the starting site then the directed path update already 
ends. The configuration is returned without being updated. Otherwise the path goes out 
through the chosen bond. As it leaves the active block it updates it by adding a monomer on 
the incoming site and then removing the monomer again if the path terminates immediately, 
or reducing the bond number on the outgoing link if the path continues. Let us invent a 
notation to describe this process. We use the symbol ® next to a bond or a site to indicate 
in-coming path into the block and use © to indicate the out-going path. For example 
(n®, bi, &_i0, bd, b_d) means that the path came into the block through the site and left 
the block through the /i = — 1 bond. When the path leaves the (active) block, the updated 
block is given by {n + 1, bi, — 1, bd, b-d)- Figure El shows a typical directed path through 
a three dimensional block. 

When the path leaves an active block through a bond it enters the neighboring passive 
block through the same bond. It is then forced to leave through one of the remaining 2d — l 
bonds with probabilities given below. Since now the value of t"^ is important, we distinguish 
two cases: < 1 and > 1. If the incoming path is through the fi = 1 (yU = —1) bond, the 
outgoing path can either be the /i = —1 (/U = 1) bond or any one of /i 7^ ±1 bonds. The 
probability to pick /i = — 1 (yU = 1) is given by 



4-2 



< 1 



2d-2+t^ 

(20) 



2t^-l 



e > 1 



K 2d- 3+2*2 

while the probability to choose any one of the fi 7^ ±1 is given by 

2d-2+fi — 



Pi,, = P-u. = ■ (21) 



L 2d- 3+2*2 t ^ J- 
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FIG. 2: A typical directed path passing through a 3d block. The path enters the block through 
the fj, = 3 bond and leaves through the /J- = I bond. The block and the path are represented by 
(n, 6i0, 62, b-2, ftsfX", 6-3)- 

However, if the path comes into the block through one of the fi 7^ ±1 bonds, then the 
outgoing path is chosen to be the /i = lor/i = — 1 bond with probabihty 



2d-2+t^ 



M-1 



< 1 



> 1 



(22) 



2d-3+2t2 

and is chosen to be one of the u 7^ ±1, bonds with probabihty 



p 



{2d-2+fi){2d~3) 



< 1 



> 1 



(23) 



2d-3+2t^ 

As the path leaves a passive block, the block itself is updated by lowering the dimer number 
on the incoming bond and raising the dimer number on the outgoing bond. We note that 
in eqs. (j2nil2ni) the probability definitions for t < 1 and t > 1 can be extended to the range 
<2d — 2 and > 1/2, respectively. This is a reflection of the fact that the requirement 
of detailed balance does not determine all the probabilities uniquely. 
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When the path leaves the passive block it enters the adjacent active block through an 
in-coming bond. In this case it can now leave through one of the remaining 2d — l bonds or 
the site. Given the in-coming path to be through the bond u, the probability to choose the 
bond fi as the out-going path is 

and the probability to choose the site as the outgoing path is 

= wk' '''' 

If the path leaves through a bond then it continues into another passive block and goes 
forward as already discussed above. However, if it exits through the site then the path 
ends. In either case the active block is again updated by increasing the dimer content of the 
incoming bond by one and decreasing the dimer content of the outgoing bond or decreasing 
the monomer content of the outgoing site by one. 

Thus, we see that the algorithm always enters and exits through the sites of active blocks. 
It increases the monomer number on the incoming site by one and decreases the monomer 
number at the exiting site by one, and thus preserving the total number of monomers. Of 
course, the site at which the directed path enters may or may not be the site in which it 
leaves, although for m = these two sites are forced to be the same since the weight of a 
configuration which contains monomers vanishes. Thus, in the chiral limit, the directed path 
is always a closed loop. An important difference between active and passive blocks is that the 
paths always enter and leave a passive block through a bond. Hence the monomer content 
of sites belonging to passive blocks never change during a single directed path update. The 
paths always enter a passive block on a bond which has at least one dimer on it. 

It remains to be shown that eqs. ()18ll2'K|l satisfy detailed balance with respect to the 
Boltzmann weights given in eqs. ()lfj|l and ()17|1 . First consider the active block represented 
by (n, 61, .., 6^, .., 6j,(8>, fe-d) with the in-coming path specified to be the bond z/. 
The probability to choose the bond ^ as the out-going bond is given by eq. p^ . After the 
update the active block is changed to (n, bi, ..,6^ — 1, .., + 1, ...bd, In order to 

prove detailed balance one is interested in the reverse process. For this, one should start 
with the configuration (n, 61, .., 6^ — l®,..,bu + 1, ...6^, where now the in-coming 
path is the bond /i, and figure out the probability of choosing u as the out-going path. This 
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turns out to be 



= .r \ , . (26) 



We see that these probabihties and the weights given in eq. (pi)|l satisfy detailed balance: 

(27) 



{N-b^)\{N-b,)\ _ b, + l {N-b^, + l)\{N-b,-l)\ 



(N-b,) b,\ bj (N-b^ + 1) ib,-l)\ (6. + 1)! 

In the above detailed balance equation we have only considered factors of the Boltzmann 
weight that change during the update. 

In the case of the passive block it is again easy to show detailed balance. For example 
consider the case when the incoming path is through the u = 1 bond and the outgoing bond 
is through the /i 7^ ±1 bond. We can represent this path by (n, 61®, 6„i, ..,b^&, ...bd,b-d)- 
After the update the new weight of the block is equal to the old weight divided by t^. 
The reverse process, starting from the updated block looks like {n,bi — l©,6_i,..,6^ + 
1®, ...bd, b-d). The detailed balance requires 

Pi.M^' = P^,i (28) 

which is indeed true when we use eqs. fl21l22j) . Using a similar method one can show that 
other choices of incoming and outgoing paths on both active and passive blocks obey detailed 
balance locally at each update. There is one exception which we discuss below. 

Consider the case when the in-coming path is a site and the out-going path is a bond and 
its reverse process. We represent this case by (n®, 61, b_i, .., 6^0, .., b^, ...bd, b_d), for which 
the probability is given by eq. ()19p . The reverse process on the other hand can be represented 
by {n + 10, 61, .., (6^ — 1)0, .., by, ...bd, b-d), which has the probability 

C = N-l, + l ^^^^ 

We now see that the detailed balance is not quite satisfied locally, since 

b, {N-b,y. N\ n + 1 (N-b, + 1)1 N\ 

N b^\ n\^{N-b^ + l) {b^-iy. (n + 1)! ^ 

There is a extra factor 1/A^ that remains uncanceled while we go from the site to the 
bond. However, since the complete "directed path" update starts and ends on an active 
site, this extra factor must appear in each direction after a complete path update. This then 
guarantees that the complete forward and reverse directed paths indeed satisfy detailed 
balance. 
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Before we end this section let us briefly comment on the ergodicity of the algorithm. 
When m — the "directed path" update is ergodic by itself. In order to see this, let us 
assume that A'^ = 1. In this case, any configuration can be transformed into any other 
configuration by a series of disconnected directed loop updates. The loops themselves can 
be identified by superimposing one configuration over the other. All dimers that differ in 
the two configurations connect the sites into disconnected loops. Clearly, there is a non-zero 
probability for performing this specific series of directed loop updates and hence to go from 
any configuration to any other. For N > 1, one can easily extend this argument and prove 
ergodicity. 

Since the directed path algorithm does not change the number of monomers, we need 
to supplement it with another algorithm that can change the number of monomers in a 
configuration when m is nonzero, in order to make the algorithm ergodic. Many options 
are available. One can choose a local algorithm based on either a Metropolis or a Heat 
Bath update. In the case of a Metropolis update for example, a simple algorithm would 
be to pick a bond at random and propose to either break the dimer into two monomers or 
create a dimer from two monomers. This proposal is then accepted with a Metropolis accept 
reject step. We find that this combined algorithm is reasonably efficient. We are currently 
exploring a more natural extension of the directed path algorithm to the massive case which 
avoids the additional local Metropohs step. This will be published elsewhere. 



IV. OBSERVABLES 



Let us for the moment assume that the directed path algorithm discussed above can 
efficiently sample configurations that contribute to the partition function. If these config- 
urations also contribute to an observable, then we may be able to compute the average 
of the observable efficiently. However, for this to be true, the observable must get all the 
contribution only from the ensemble of configurations sampled by the algorithm and the 
value of the observable must not fluctuate much in this ensemble. Observables which satisfy 
these properties will be referred to as "diagonal" observables. For example observables like 
dimer-dimer correlations and spatial winding numbers associated to the global U{1) chi- 
ral symmetry are diagonal observables. When m is not very small, monomer correlations 
can also be treated as diagonal observables. On the other hand when m — 0, observables 
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involving monomers do not fall in this category; such observables get contributions from 
configurations with monomers, while the algorithm only produces zero-monomer configura- 
tions. In other words the configurations that are important to the partition function are not 
useful to measure observables. Such observables are difficult to compute and will be referred 
to as "non-diagonal" observables. The chiral condensate and the chiral susceptibility are 
two interesting examples of such observables. We will argue below that the directed path 
algorithm offers an efficient method to compute them. 

From the discussion in the previous section we know that a directed path update starts 
on an active site, goes through an alternate series of both passive and active sites and ends 
on an active site. It is important to recognize that between the start and the end of each 
update we sample a configuration with exactly two additional monomers every time we 
visit a passive site; one monomer is located at the starting active site and the other at the 
visited passive site. These intermediate configurations can be used to compute monomer 
correlations. To write down explicit relations we introduce an integer function I{x,y) for 
each directed path update where x, y are lattice sites. This function is defined as follows: 
Before we start the directed path update we set I{x, y) = for all values of x, y; we then 
add one to I{x,y) if the directed path update starts on the active site y and every time it 
visits the passive site x. Now consider a monomer-dimer configuration with Ui monomers 
located at i = 1,2, k. For this configuration we can define a one point function and a 
two point function as follows: 

^'(^) = 2l^J^' S /(-.^) (31) 



S2{x,y) ^ 



\i=l 



(32) 



2d -2 + 2^2 
It is possible to show that 

i( ) = (Srix) ) (33) 

^{^i;{x)^^{y)) = {S2{x,y)). (34) 
In appendix |X1 we give a brief derivation of these results. Using eqs. and one finds 



{^^^) = {J2Sl{x)) (35) 

X 

and 

X = ({n + l)J2S2ix,y)), (36) 
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where we have used n = ni + n2 + ■■■ + Uk to denote the total number of monomers in the 
configuration. In eqs. ()33II36|) the average on the right hand side is taken over the ensemble 
of configurations generated in the directed path algorithm. 

As the derivation in the appendix shows, some additional work is necessary to derive 
expressions for the non-diagonal observables in the directed path algorithm. However, it is 
usually possible to find expressions for most observables. In fact similar expressions have 
been obtained in the context of other cluster algorithms 0,^^. In the next section we will 
demonstrate the correctness of the relations (jH^j) and (jHUj) by comparing the results obtained 
using them with exact calculations to a high accuracy. 

It is clear that one can also compute many other interesting quantities. All observables 
involving two monomers can be obtained using eq. (j34|) . Higher point monomer correlations 
functions can also be computed but need additional work. They will be discussed elsewhere. 



V. TESTING THE ALGORITHM 



We have tested the directed path algorithm and the expressions of the chiral condensate 
fl35|) and the susceptibility (jHUj) in various dimensions for small lattices where exact calcu- 
lations are possible. In this section we briefiy review our tests. For a finite lattice with V 
sites, the partition function (eq.®) is an even polynomial in the mass 

Z{m) = Co + C2m^ -|- C4m^ + . . . + CNvTri^'^ ■ (37) 

The coefficients C2n are functions of t^, although we suppress this dependence in the notation. 
The absence of terms with odd powers of m in eq. ()37|) is a consequence of the remnant U{1) 
chiral symmetry. Every configuration [n, b] can only have an even number of monomers. 
NV is the maximum number of monomers allowed. 

The condensate and the susceptibility can be obtained from the partition function using 
the relations 

- 1 1 dZ{m) _ 1 1 d'^Z{m) 

~ V ^{m) ~d^ ' Y^m) dm? 

Thus, once the coefficients cq, C2, . . . , cnv are known these quantities can be determined. In 
the appendix we give expressions for these coefficients in some simple cases. When m = the 
condensate vanishes as it should and x = 2c2/l^Co. Hence, in some cases where calculating 
all the coefficients is difficult, we give only the coefficients Cq and C2. 
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In tables m and ini we compare the condensate and susceptibihty obtained using and 

with exact results. 

TABLE L Chiral Condensate: Algorithm vs. Exact results. 

chiral condensate 



N 


Lattice Size 


m 




Exact 


Algorithm 


1 


2x2 


0.3 


1.3 


0.12133... 


0.12135(4) 


1 


2x2x2 


0.47 


0.43 


0.716859... 


0.71658(23) 


2 


2x2x2 


0.1 


3.6 


0.067350... 


0.06736(2) 


1 


2x2x4 


0.2 


6.4 


0.0206279... 


0.02063(3) 


1 


4x2x2 


0.4 


4.7 


0.1263817... 


0.12643(5) 



TABLE II: Chiral Susceptibility: Algorithm vs. Exact results, 
chiral susceptibility 



N 


Lattice Size 


m 




Exact 


Algorithm 


1 


8 X 


8 


0.0 


1.0 


5.27221... 


5.2722(2) 


3 


4 X 


4 


0.0 


1.0 


14.1595... 


14.159(8) 


30 


2 X 


2 


0.0 


1.0 


338.534... 


338.2(8) 


1 


2x2 


X 2 


0.13 


1.2 


0.57028... 


0.5702(2) 


3 


2x2 


X 2 


0.0 


1.4 


4.38869... 


4.3884(13) 


1 


2x2 


X 4 


0.0 


3.2 


0.25766... 


0.2576(1) 


2 


2x2 


X 4 


0.0 


1.4 


3.17378... 


3.1741(9) 


1 


4x2 


X 2 


0.0 


0.37 


0.69463... 


0.6945(3) 


2 


4x2 


X 2 


0.0 


3.2 


2.24205... 


2.2418(5). 


1 


2 X 2 X 


2x2 


0.0 


5.3 


0.15401... 


0.15397(6) 


1 


2 X 2 X 


2x2 


0.0 


1.2 


0.80231... 


0.8022(2) 



In order to check if the algorithm performs well for large lattices, we have derived an 
exact expression for x on a 2 x L lattice when m = 0, = 1, t = 1. After a lengthy 
calculation we found 

^ ^((V2 + ir^ + (V2-ir^)-l 
^ (V2 + l)^ + (72 -1)^ + 2 ^ ^ 
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For L ^ oo this becomes 

lim X = ^^^^^ = 0.8535534... (40) 
2v2 

Our algorithm yields x = 0.8535(3) when L = 1024. 

VI. PERFORMANCE OF THE ALGORITHM 

An important feature of a good algorithm is a short auto-correlation time which we 
denote as r. Due to critical slowing down the auto-correlation time increases with increasing 
correlation lengths. One expects r (x where C, is the correlation length and z is the 
dynamical critical exponent of the algorithm. For most local algorithms 1 < z < 2. Many 
efficient cluster algorithms on the other hand are known to have < z < 1. In this section 
we estimate z for the = 1 algorithm in two dimensions for m = and = 1. 

Although the auto-correlation time depends only on the algorithm, a more useful quantity 
in practice is the integrated auto-correlation time Tint defined for a given observable using 
the relation 

r^nt = - E G{t), (41) 



2 



where 



with Oi being the value of the observable generated by the algorithm in the i^^ sweep and (O) 
refers to its average. Ideally when the correlation function G{t) is known exactly then tmax 
can be set to infinity. However, in practice due to a finite data sample, G{t) is determined 
with growing relative errors as t increases. Hence, tmax must be chosen carefully jl3|. It 
is also important to normalize Tint such that the definition of a sweep does not enter into 
it. Here, every directed loop update after a site is picked at random is defined as a sweep. 
This means that on an average it takes as many as / = L^/ (A/'passive) sweeps ( A/'passive is the 
number of passive sites encountered in the directed path update) before we update every 
degree of freedom once. Hence we divide the answer obtained from eq. by / and define 
that as the integrated auto-correlation time. 

Here we estimate Tint for the chiral susceptibility. As we will see in the next section, 
^ is infinite for the parameters we use (A^ = l,d = 2,m = 0), since there are massless 
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FIG. 3: The auto- correlation function G{t) is shown as a function of t for L = 16, 24, 32, 48 and 
64 

excitations. Hence, we study the behavior of 

^int ^ 

function of the lattice size L. Figure 
El plots the function G{t) for L = 16,24,32,48 and 64. This graph suggests a simple choice 
for tmax for various L's. We choose tmax = 12, 16, 20, 30, 30 for the five values of increasing 
L. This defines Tint uniquely. We plot the results in figure El The function 0.107L°-^'^ (solid 
line) roughly captures the behavior of function of L suggesting that the dynamical 

critical exponent z of our algorithm is around 0.5. 

Let us now compare the directed-path algorithm discussed here with the Metropolis 
algorithm developed in a previous work 10|, lll|. This is of interest since the results from 
the previous algorithm were quite unexpected and may be wrong. In the earlier work the 
two dimensional partition function was rewritten in terms of loop variables and the dimers 
along the loop were updated using a Metropolis accept-reject step. It was shown that the 
algorithm reproduced exact results with good precision on small lattices. Based on the 
evidence from the directed path algorithm, we now argue that the auto-correlation time 
of the previous algorithm grows uncontrollably for large lattices. In figure El we compare 
4000 consecutive measurements of the chiral susceptibility in the simulation time history 
at m = on a 32 X 32 lattice, between the old and the new algorithms. In both cases. 
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% 10 20 30 40 50 60 70 

L 

FIG. 4: Integrated auto- correlation time as a function of L. 

measurements are performed after a comparable amount of computer time. The fluctuations 
in the old algorithm are clearly much larger. The average of 300, 000 measurements with 
the old algorithm gives a value of 41.8(6) as compared to 43.669(5) obtained from just 
4000 measurements with the new algorithm. This clearly demonstrates the efficiency of 
the new algorithm in addition to showing that the earlier algorithm underestimates the 
final answer. Note that there is one measurement of the order of 10^ visible in the old 
algorithm time history. If we remove this measurement from our analysis we obtain an 
answer of 41.0(2) even with 4000 measurements. This is consistent with the value obtained 
by averaging over the entire time history. Without more statistics it would be impossible 
to say whether these spikes are isolated events or contribute an important fraction to the 
final answer. However, our new and more efficient algorithm shows us that the contributions 
from these isolated spikes are indeed important and can contribute up to five percent to the 
final answer on a 32 x 32 lattice. Comparing the results obtained from the new and the 
old algorithms for different lattice sizes, we now believe that these large spikes contribute 
even a larger percentage to the final answer on larger lattices. Most likely the old algorithm 
is not able to tunnel between important sectors of the configuration space and develops 
large auto-correlation times at large volumes. Due to insufficient statistics we think that 
the contributions from the spikes were not sampled properly in the previous work and hence 
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4000 



%44 




4000 



FIG. 5: Comparison of simulation time histories of the chiral susceptibility for 4000 consecutive 
measurements between the old (top graph) and the new (bottom graph) algorithms for L = 32. 

the final answers were systematically lower. This led to wrong conclusions. 



VII. RESULTS IN THE CHIRAL LIMIT 



In this section we present some results obtained using the algorithm in the chiral limit. We 
focus on the issue of chiral symmetry breaking in various dimensions. As discussed earlier, 
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the finite size scaling of tlie cliiral susceptibility is an ideal observable for investigating this 
issue. All the results given below are obtained with t = 1. 



le+06 
le+05 
10000 

X 

1000 
100 



8 10 12 16 24 32 48 64 96 128 192 

L 

FIG. 6: Finite size scaling of the Chiral Susceptibility in three and four dimensions. The solid 
lines represent fits of the data as described in the text. 

Based on the work of we know that at strong couplings chiral symmetry is broken in 
three and four dimensions. In figure IHl we plot the chiral susceptibility (%) as a function of 
the lattice size L in ci = 3, 4 when = 1 and in ci = 4 when = 3. For = 4 the data fits 
very well to the functional form aL'' + 6, where a and h are constants. This functional form 
is motivated from chiral perturbation theory. In the case of = 3 chiral perturbation theory 
suggests a form aL'^ + bL'^~'^ + c. Again the data fits well to this form over the entire range. 
The values of all the fit parameters and the x^/d.o.f. obtained from the fits are given in 
table ITTTl The solid lines in figure El represent these fits. The divergence of the susceptibility 



TABLE III: 



N 


d 


a 


b 


c 


xVd.o.f. 


1 


3 


0.05493(5) 


1.75(9) 


-10(1) 


0.5 


1 


4 


0.05124(7) 


3.6(2) 




0.5 


3 


4 


0.4823(2) 


12(5) 




0.2 
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with L is consistent with spontaneous breaking of chiral symmetry. Using eq.()14|) we see 
that the infinite volume chiral condensate is given by = \fa. We find 



0.2343(2) for N =\,d = Z 
0.2264(2) for N =\,d = ^ 
0.6945(3) for = 3, = 4. 



(43) 



We would like to remind the reader that the errors in the fitting parameters do not include 
all the systematic errors that arise due to variations in the analysis procedures. 

In two dimensions the Mermin- Wagner- Coleman theorem forbids the formation of a con- 
densate jisl . Q|. However, it is well known that a theory with f/(l) symmetry can still 
contain massless excitations. Of course there is always a possibility for the theory to be in 
the chirally symmetric massive phase. To which phase does our model at t = 1 belong? The 
behavior of x as a function of L in both the possible phases was already discussed in eq. lfT^ . 
In figure [7| we plot % as a function of L for = 1, 2, 3, 5, 8, 10, 20. We find that all of the 
shown data fit reasonably well to the function x = aL?'^^ . The solid lines show these fits. 
The values of a, rj and x^/d-o.f. are given in table II VI This result strongly suggests that two 
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FIG. 7: Finite size scaling of the Chiral Susceptibility in two dimensions for different values of N. 
The solid lines are fits as discussed in the text. 

dimensional staggered fermions at strong couplings are in the critical massless phase. 
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TABLE IV: 



N 


a 


7] 


xVd.o.f. 


1 


0.239(1) 


0.498(2) 


1.2 


2 


0.586(1) 


0.230(1) 


1.6 


3 


1.139(3) 


0.149(1) 


2.7 


5 


2.78(1) 


0.086(1) 


0.11 


8 


6.71(2) 


0.054(1) 


0.28 


10 


10.22(5) 


0.043(2) 


0.33 


20 


39.1(2) 


0.021(2) 


0.64 



The partition function of two dimensional dimer models (A^ = 1, no monomers) on a 
planar lattice can be solved exactly [2^ 13, • In [2^ a determinant formula for the two 
monomer correlation function ((■?/''?/' (x) ipip{y))) was derived in the infinite lattice setting. A 
combination of analytic and numerical analysis in 2^ provided strong evidence that this 
correlation function decays as l/|x — ?/|^/^ when |x — ?/| is large. This implies (in = 1 case) 
that the chiral susceptibility x diverges as L^/^ in large volumes, i.e. rj = 0.5. Our N = 1 
algorithm result 7] = 0.498(2) is in excellent agreement with this semi-exact result. 

At the other extreme, a large N analysis leads to the conclusion that chiral symmetry is 
spontaneously broken even in two dimensions. Then the condensate must be non-zero. How 
does the condensate acquire a non-zero value? A resolution to this paradox was proposed 
by Witten in He predicted that rj c/N for large A^, so that when N is strictly infinite 
7] = and X ~ as expected in a phase with broken chiral symmetry in two dimensions 
(see eq. fll4j) ). Our data is consistent with this expectation. We find that all of our data fits 
well to the form r] = 0.420(3)/A^ + 0.078(4) /A^^ ^j^j^ xVd-o.f. ~ 0.33. Our data along with 
the fit is shown in figure |H1 The 1/A^^ term in the fitting function is necessary to fit the data 
for A^ = 1,2,3. 



VIII. CONCLUSIONS 



In this article we have introduced a new type of cluster algorithm to update strongly 
coupled U{N) lattice gauge theories with staggered fermions in the chiral limit. We studied 
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FIG. 8: 7] as a function of N. The solid line corresponds to the function 0.42/iV + 

the performance of the algorithm in two dimensions for a specific set of parameters and 
discovered that it has a dynamical critical exponent of roughly one half. We found the 
algorithm to be extremely efficient even in three and four dimensions. This progress allows 
us, for the first time, to study questions related to the chiral limit of a gauge theory beyond 
the mean field approximation. 

As a first step we studied the question of chiral symmetry breaking. We used finite size 
scaling of the chiral susceptibility as a tool to address this question. We found clear evidence 
that the U{1) chiral symmetry is spontaneously broken in three and four dimensions. On 
the other hand, in two dimensions we found that the theory contains massless excitations al- 
though the chiral condensate was zero. This is consistent with the Mermin-Wagner-Coleman 
theorem. The critical exponent 1] = 0.498(2) obtained with the algorithm for = 1 in two 
dimensions is in excellent agreement with the semi-exact result of 0.5 inferred from 1251 1. 

n 

Further, r] scales as for large as predicted in |21l |. 

We believe our algorithm can be useful in addressing some interesting physics questions. 
For example, the universality of the finite ternperature QCD chiral phase transition with 
staggered fermions remains controversial Q|- It would be useful to give a definitive 
answer at least at strong couplings. Previous calculations could only be performed far from 



the chiral limit and on small lattices 



12| |. Using our algorithm we can now settle this question 
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completely. Another question of interest is to determine the critical density of baryons where 
chiral symmetry restoration takes place. At strong couplings and in the static approximation 
(heavy baryon limit) this question can be answered with our algorithm. Since baryons are 
known to be heavy at strong couplings this approximation may even be justifiable. Our 
algorithm can also be used to extract quantities like the decay width of the rho meson at 
strong couplings. This would be instructive since we would learn more about the difficulties 
involved in such calculations that are not merely algorithmic. 

The present work also raises many interesting algorithmic questions. Can one extend 
our algorithm to other situations? For example, can one study more than one flavor of 
staggered fermions? This is interesting since one can then explore more complex chiral 
symmetry breaking patterns. Can Wilson fermions and Domain wall fermions be studied 
at strong couplings using similar techiiiques? Interesting and controversial phase structures 
have been predicted for the latter S^lsOl- What about a systematic way to go towards the 
weak coupling limit? Perhaps allowing a fermionic plaquette term in the Boltzmann weight 
already has some desirable effects. What about more general applications? For example it is 
possible to rewrite the determinant of a matrix as a sum over monomer-dimer configurations. 
Does this mean we can calculate the determinant of certain types of matrices more easily 
than before? 

Finally, our algorithm should also be of interest to physicists studying monomer-dimer 
systems. Such systems are interesting in their own right and have a long history [SI]. Many 
statistical mechanics problems including the Ising model can be reformulated in terms of 
monomer-dimer models on various kinds of lattices. Novel Monte Carlo algorithms continue 
to be developed to study them js^. We believe that our approach can provide a useful 
alternative. 
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APPENDIX A: DERIVATION OF EXPRESSIONS FOR OBSERVABLES 



Here we present a brief derivation of the relations ()33p and ()34|) . We begin by noting that 
the chiral condensate is given by 

(^^(x)) = -l- W{[n,bl)m- (Al) 

where the sum is over monomer- dimer configurations [n, h]x which satisfy the usual constraint 
relation on all sites except the site x where the relation is modified to 

nx + Y.^^,{x)+h^^{x) = N -I. (A2) 

The partition function on the other hand is given by 

Z{m) = W^(K6].)m" (A3) 

[n,b] 

where the sum is over monomer-dimer configurations [n, h] which satisfy the usual constraint 
on all sites. The weights iy([n, 6]^) and appearing in ()A1|) and ()A3|) are obtained 

by multiplying the block weights ()16|) and ()17|) associated to all active and passive blocks. 
Since the mass is not taken into account in these weights we have an extra factor where 
n refers to the total number of monomers in the configuration. 

Let us now construct an update to go from a given configuration [n', h']^ to a configuration 
[n, h]. This update is performed as follows: 

1. We declare the site x to be passive. 

2. We start at the site x and choose the first bond to be one of the 2d bonds. The 
probability to choose the /i = ±1 bond is given by t^/{2d-2 + 2t^) and the probability 
to choose one of the 7^ ±1 bond by l/{2d — 2 + 2t^). 

3. Once we pick the bond we increase the dimer content of that link by one and go to 
the adjoint active block. 

4. We then use the probabilities discussed in section IIIII to construct a directed path 
update which ends on some active site y. 

Let the path generated using the above rules be referred to as a^y At the end of the 
above update it is easy to prove that the new configuration belongs to the type Let 
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P{\n' , h']x] axy) refer to the probability to produce the path axy starting from the configura- 
tion [n' , h']x using this update. 

Now consider the unique reverse path of a^y referred to as a~y. This path is one of the 
"partial" directed paths that we produce using the directed path algorithm described in 
section 11111 In particular this path is a path that starts from the active site y and visits 
the passive block at x. Let P{\n,h]]a~y) be the probability of generating this reverse path 
using the algorithm. Since the forward and reverse paths are produced by probabilities that 
satisfy detailed balance at each stage except for possible factors at the ends, it is easy to 
argue that 

P(K, h%- axy) W{[n\ h']x) m"' = _ ^ ^ P([n, 6]; <) W{[n, h]) (A4) 

The factor V arises because the site y is picked with probability 1/V but not the site x. The 
factor is due to the uncanceled factor from (jHIH) . The factor l/{2d — 2 + 2t^) compensates 
the same factor from the left hand side. The mass factor arises because of the mismatch in 
the number of monomers between the two configurations; in particular n' = n + 1. 
Now by construction we know 

^ P{[n\h']x-axy) = l (A5) 

where the sum is over all possible paths axy These paths always start at the same x but 
end at various sites y. Thus we derive the relation 

where the configurations [n, b] on the right hand are determined from the configuration 
[n', b']x and the path axy It is then possible to argue that 

E W{[n',b']x)m-' = E E (2d-2 + 2t^) Pi[n,b];ay^')Wi[n,b])m\ (A7) 

Since the directed path update produces paths a~y starting from the configuration [n, b] 
with probability P{[n, b]; a~y), we see that 
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where I{x,y) was defined in section HVl and the average on the right hand side is taken 
over the ensemble of configurations generated in the directed path algorithm. This proves 
relation 

In order to show (j34p we start with the relation ()A6|) and assume that the mass factors 
are site dependent. By differentiating the relation ()A6|) with respect to the mass factor at 
the site y we can generate weights of configurations that contribute to the two monomer 
correlations. This then yields ()34|) . We leave the steps of this derivation to the reader. 

APPENDIX B: EXACT RESULTS ON SMALL LATTICES 

In this appendix we give explicit expressions for the coefficients C2n defined in eq. ()38p. 
as a function of x = t'^ for a selection of small lattice sizes and small values of where exact 
calculations are possible. 

1. iV = lona2x2 lattice 

Co = 4(1 + 
C2 = 4(1 + a;) 

C4 = 1 

2. iV = lona2x2x2 lattice 

Co = 16(4 + 4x2 + 

C2 = 8(16 + 16a; + 8a;2 + 4x^) 

C4 = 4(20 + 16a; + 6x'^) 

C6 = 2(8 + 4x) 

Cg = 1 
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3. A/^ = lona4x4 lattice 

Co = 16(1 + Ax'^ + 7x^ + 4:x^ + x^) 

C2 = 64(2 + 4x + lOx^ + 13x^ + 13x^ + IQx^ + 4x^ + 2x^) 

C4 = 32(13 + 40x + Slx^ + 96x^ + 81x^ + 40a;^ + 13x^) 

C6 = 64(11 + 37a; + 63a;2 + 37x^ + llx^) 

cg = 8(83 + 256x + 354^2 + 256x^ + 83x^) 

Cio = 32(11 + 28x + 28x2 + llx^) 

ci2 = 8(13 + 24a; + 13x2) 

Ci4 = 16(1 + a;) 

Cl6 = 1 

4. iV = 2ona2x2x2 lattice 

Co = 1156 + 3136a;2 + 2116x^ + 576a;^ + 81a;^ 

C2 = 16(476 + 784a; + 1064a;2 + 870a;^ + 502a;^ + 264a;^ + 72a;^ + 27a;^) 
C4 = 4(4428 + 8512a; + 9532^2 + 6624,t-^ + 3302x^ + 1056a;^ + 243a;'') 
C6 = 16(1132 + 2116a; + 1976^2 + 1116a;^ + 386a;^ + 75x^) 
Cg = 2(4714 + 7744a; + 5768a;2 + 2304a;^ + 443a;^) 

Cio = 16(166 + 220a; + 116a;2 + 25a;3) 

ci2 = 12(34 + 32x + 9x2) 

ci4 = 16(2 + x) 

Cl6 = 1 

5. A/^ = lona4x2x2 lattice 

Co = 16(256 + 256x2 + 112x^ + 16x^ + x^) 
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C2 = 128(128 + 128x + 160x^ + lOAx^ + 52x^ + 20x^ + Ax^ + x'^) 

C4 = 32(832 + 1280x + 1296a;^ + 76Sx^ + 324x^ + 80x^ + 13x^) 

ce = 64(352 + 592X + 504x2 + 252x3 + 74a;^ + llx^) 

cg = 8(1328 + 2048x + 1416x2 + 512^3 + 83x^) 

cio = 32(88 + 112X + 56x2 + 11x3) 

ci2 = 8(52 + 48x + 13x2) 

Ci4 = 16(2 + x) 

Cl6 = 1 

lona2x2x4 lattice 

Co = 16(81 + 232x2 + 216x^ + 96x'^ + 16x^) 

C2 = 128(63 + 98x + 142x2 + 130x3 + 84x^ + 56x^ + 16x^ + 8x^) 

C4 = 32(563 + 1064x + 1242x2 + 944x3 + 532x^ + 192x^ + 56x^) 

C6 = 64(284 + 529x + 512x2 ^ 3^2x3 + 120x^ + 28x^) 

cg = 8(1179 + 1936x + 1492x2 + 640x3 + 140x^) 

cio = 32(83 + 110X + 60x2 + 14x3) 
= 8(51 + 48x + 14x2) 

Ci4 = 16(2 + x) 

Cl6 = 1 

2 on a 4 X 4 lattice 

Co = 65536(81 + 576x2 ^ 2416x^ + 5648x^ + 7520x^ + 5648x^° 

+2416x^2 ^ 57g2,i4 ^ g^^ie^) 

C2 = 2097152(54 + 144x + 564x2 + 1145x3 + 2490x^ + 3806x^ 
+5470x^ + 6303x^ + 6303x^ + 5470x^ + 3806x^° + 2490x^^ 
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8. Af = 3ona2x2x2 lattice 

Co = ^(4624 + 25600x^ + 35396x^ + 20224x^ + 5800x^ + 900x^° + 81x^2) 
C2 = ^(10880 + 25600a; + 54080a;2 + 63556a;^ + 629120;^ + 47560a;^ 



9. A/^ = 2ona4x2x2 lattice 

Co = 1336336 + 3625216x^ + 5242816^^ + 3817024x^ + 1491668x^ 

+318272x^° + 37072x^2 ^ 2304^1^ + 81x^^ 
C2 = 64(275128 + 453152x + 1053864^^ + 1260996x^ + 1522636a;^ 
+1299408^^ + 1009344x^ + 630309x^ + 334475x^ + 154496^^^ 

+56066a;^° + 19137a;^^ + 4508a;^2 + 1128a;^^ 



10. AT = 2 on a 2 X 2 X 4 lattice 

Co = 198916 + 1682272^^ + 4177396^^ + 4825728^^ 

+3184704a;^ + 1343232a;^° + 381996x^^ + 69984a;^^ + 6561a;^^ 

C2 = 32(140936 + 399424a; + 115119682;^ + 1787069^;^ + 2560377^;^ 
+2707684a;^ + 2511784a;^ + 2002594a;^ + 1320438^^ + 840288a;^ 
+403668x^° + 212841a;^^ + 70956x^^ + 31590a;^^ + 5832a;^^ 



+28220a;^ + 15120x^ + 5620a;^ + 2220x^ + 450a;^° + 135a;^^) 



(Bl) 



+144x^^ + 27a;^^) 



(B2) 



+2187a;^^) 



(B3) 
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11. Ar = lona2x2x2x2 lattice 



Co = 256(81 + 124x2 + 54x^ + 12x^ + x^) 
C2 = 128(792 + 968x + 936x2 + 600x^ + 240x^ 
+144x^ + 24x^ + Sx'^) 
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